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文章 编号 :1005-3085(2010)04-0652-11 


形状 梯度 法 求解 不 可 压缩 流 的 形状 优化 问题 * 


朱 静 涛 ， 高 志明 ， HRE 
(西安 交通 大 学 理学 院 ， 西 安 710049) 


摘 要 : 流体 中 物体 形状 优化 设计 在 实践 中 有 重要 的 应 用 。 对 于 区 域内 由 Navier-Stokes 方 程 描述 的 流 
体 ， 本 文 研究 以 流体 状态 的 泛 涵 为 目标 函数 的 优化 问题 。 基 于 共 轿 方法 与 函数 空间 参数 化 方法 ， 
本 文 得 到 了 问题 的 形状 导数 。 在 此 基础 上 构造 了 一 种 共 轿 梯度 算法 。 数 值 例 子 表 明 本 文 的 方法 是 


可 行 的 和 稳定 的 。 
关键 词 : EARR: Navier-Stokes 方程 ， 极 大 极 小 值 问 题 
分 类 号 : AMS(2000) 65K10; 76B55; 49Q20 中 图 分 类 号 : 0241.82; 0242.21 文献 标识 码 : A 
1 引言 


不 可 压缩 流 的 形状 优化 问题 在 很 多 领域 都 有 着 广泛 的 应 用 ， 如 航空 、 发 动机 推进 、 水 利 、 
海洋 学 、 建 筑 学 、 空 气动 力学 等 中 。 本 文 的 目标 是 根据 Lagrange 方 程 的 极 大 极 小 值 问题 的 可 
微 性 来 推导 关于 目标 函数 的 共 斩 方程 。 再 使 用 函数 空间 参数 化 方法 推导 关于 目标 函数 的 形状 导 
数 。 使 用 形状 导数 来 求解 二 维 混合 边界 条 件 的 不 可 压缩 流 的 形状 优化 问题 。 

本 文 的 结构 内 容 如 下 。 第 二 部 分 给 出 不 可 压缩 流 的 形状 优化 模型 和 速度 法 介绍 。 第 三 部 分 
推导 对 于 给 定 的 目标 函数 的 共 斩 方 程 。 第 四 部 分 使 用 函数 空间 参数 化 方法 推导 出 关于 给 定 目标 
函数 的 形状 导数 公式 。 第 五 和 六 部 分 分 别 给 出 了 形状 优化 的 具体 算法 和 数值 算 例 ， 验 证 了 算法 
的 有 效 性 和 稳定 性 。 


2 ”预备 知识 


首先 给 出 不 可 压缩 流 的 形状 优化 问题 模型 。 设 是 RN (N = 2,3) 上 的 有 界 开 区 域 ,，Q C 
日 是 开 区 域 ( 见 图 1D)， 它 的 边界 记 为 839， 它 是 C2 类 的 。 如 图 1，Fin 是 入 流 边 界 ，Tout 是 出 流 
WH, To = 6990/(Tin UTout)。 从 Tin 给 定 一 个 入 流 g， 经 过 管道 Q 最 后 从 Tout 流 出去。 在 管 
壁 To 上 流速 度 以 为 0。 


—vAut+(Du)-u+Vp=f, E Q P, 


divu = 0, 在 Q 中 ， 

u=g, 在 _ Pin E, (1) 
u=0, 在 To E, 

o-n=0, 在 Fout E, 
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其 中 = 1/Re，Re 为 Reynolds 数 ,0o = —pl + v(Du + D*u), u = (wi,w2); 


Du = ( Oyu, iuz ) , (2) 
zuı zuz 
D*w 为 Du 的 转 置 。 
令 Tp := Tin To, WutElp kA u=Qg, 其 中 
》 在 Tin E 
g=4 7 (3) 
0, #2 To EE. 





图 1: ARF RRQ 


对 于 给 定 优 化 问题 的 目标 函数 .J(Q)， 我 们 求解 形状 优化 问题 就 是 寻找 一 个 区 域 0* c OK 
得 
J(9*) = min J(Q):= min F(Q, u, p), (4) 
E (u, p) WH (1) R- 
本 文中 考虑 以 下 两 种 目标 函数 的 形状 优化 问题 ， 其 中 wa 是 事先 给 定 的 已 知 的 目标 速度 ， 在 
实际 应 用 中 一 般 由 设计 者 给 出 


719) =3 f lu- valde, (5) 
HQ) = 5 [ w-uaPas, (6) 


Hp c DD 是 需要 优化 的 边界 。 以 下 我 们 主要 以 第 一 个 目标 函数 (5) 来 讨论 形状 优化 问题 的 算 
法 推导 ， 目 标 函数 (6) 讨论 类 似 。 并 在 第 六 节 中 给 出 了 对 于 这 两 类 目标 函数 的 数值 算 例 。 

这 里 先 简要 介绍 一 下 速度 法 的 理论 5656。 设 区 域 的 边界 9 是 C 的 , 上 > 1。 令 V(t,x) E€ 
EF := C((0,7]; (RY, RY)), 其 中 ORY, RNY) 是 k 阶 连续 的 具有 紧 文集 的 函数 空间 , RY 是 入 
维 欧 氏 空间 。 考 虑 以 下 动力 系统 


向 (t,X)= V(t, z(t)), t>0, 
x(0, X) = X, 


(7) 
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其 中 和 是 QQ 中 的 点 。 ST, : R" >R”, Q(X) = r(t, X) KP r(t, X) Æ t HŽ (7) RR. 
定义 


OQ := TQ) = {2 ER” |3X EQ, 使 得 z(0) = X, 2(t) = 2}. (8) 
定义 16 目标 函数 .J :0 >R 
(i) 如 果 极 限 


lim J(%) ~ TO) _ 
t 一 0 t 
存在 ， 则 称 了 在 区 域 人 上 沿 着 V 方向 上 的 欧 拉 半 导数 为 dJ(Q;V)。 

(ii) 如 果 dJ(Q;V) 存在， 对 于 任意 V c EN 映射 一 dJ(Q;V) : EN 一 及 是 线性 且 连 续 
的 ， 则 J 了 在 Qf 上 是 形状 可 微 的 ， 且 dJ(Q;V) 定义 为 J 在 Q 上 的 形状 导数 。 

由 Hadamard-Zolésio 结构 定理 [四 可 知 关 于 函数 的 Euler 半 导 数 总 可 以 化 为 以 下 形式 


:dJ(Q;V) 


dJ(Q;V) = (K(8Q), yea(V) n), (9) 
Etn EUROON LINEE, yoa BV ZEON 上 的 迹 。 由 (9) RF AAT 的 欧 拉 半 导 数 可 以 
限制 在 9Q LIF ARMS V E 上 的 限制 有 关 ， 则 9 吕 上 的 形状 梯度 为 入 (60) := yho(Kn), H 
中 yo 为 Yan 在 69 上 的 伴随 算 子 。 
3 RAE 


下 面 讨论 目标 函数 (5) 式 对 应 优化 问题 中 的 共 斩 方程 。 首 先 给 出 方程 (1) 式 的 变 分 形式 。 
设 Q 具有 光滑 的 边界 ， 且 简 记 Q 上 的 一 些 空间 如 下 


V5(2) := { (v, 1) |v € H?(Q), le H1(Q)/R, v= 9 Tp E}, 
Vo(2) := { (v, 1) |v € H?(Q), le H*(Q)/R, v =0 在 Tp E}, 


其 中 型 4(D) = (H4(D))”. 假设 f 和 g 是 固定 的 ，(w,p) 是 (1) 的 解 ，f, 9g € HRY), (up) € 
Vs(0). 
令 
alp, 4%) := fi 2ve(p) : e(W)dz, 
2 


其 中 e(p) = iDo + D*p)， 符 号 4 : BRAM AM BM NTR AM, MË A = 
(aij)mxm, B= (bij)jmxm WA: B= Dy; aijbigo 定义 
blg, 7) := f naiveds, cp; := | (Dop: dae, (p,p) := f o- pàs, 
WW (u,p) 是 (1H) 式 的 解 当 且 仅 当 满足 下 式 ， 也 就 是 N-S 方 程 的 变 分 形式 
E(Q;u,p;9,9) =0, VY (p,n) € Vo(Q), (10) 


其 中 
E(Q; u, p; p,n) := alu, p) — b(y, p) — b(u, n) + c(u;u, p) — (f, p). (11) 
这 样 (5) 和 (10) 构成 一 个 带 等 式 约 束 的 优化 问题 。 我 们 需要 在 区 域 Q0 上 通过 djJ(9;Y) 指示 
的 方向 V 变 化， 寻找 J(Q) 的 欧 拉 半 导 权 。 
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首先 ， 我 们 构造 如 下 Lagrange 函数 
G(O; p, my, 7) := J(Q) + EQ; p,n V, 1), (12) 


则 在 确定 的 区 域 Q 时 ， 只 需求 以 下 极 大 极 小 值 问题 


Boe. Ce G(Q; 9,05; p, y) (13) 
G 的 鞍点 应 满足 下 列 条 件 
Zin. b+) VEVD 04 
a a 
Ba Oi Bs Ps n)du + apo u, p;p, n)ôp=0, Y (u,p) € Vs(0). (15) 


公式 (14) 转化 为 对 应 的 偏 微分 方程 也 就 是 Navier-Stokes 方 程 (1)。 公 式 (15) 也 可 以 有 对 应 
的 偏 微分 方程 形式 。 推 导 如 下 


0 
ed u, p; p,n)ôp = —b(¢, dp), (16) 
0 
— G(Q; u, p; p, n)du = J (u — ua) - dudz + a(du, p) 
ðu Q 
— b(5u, n) + c(du; u, p) + e(u; du, p). (17) 


因为 divu = ul 十 zuz = 0, 则 


2 
c(u, du, p) = f (Diuyu- par = D | vindibujde 
2 Q 


i,j=1 


2 
= u:n)(du-y)ds 一 /ww dz+ f :6u;Oiuid 
[emu pas- D> (fswsorwsdr | 6u0mmdr) 


ij=1 


2 2 
= f „© môu: p)ds — I (Dy)u- dud — X- ( 人 pjðuj 》 druid) 


Æ E 
z [ (uu pds — elu, 9, 5u); 
Siac A ORE a (D*u)p - dude. 
那么 (15) 式 也 可 以 写成 
= f 人 


+ 人 -n)(du- p)ds — c(u; p, du) + [wwe - dudz, (18) 
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则 (18) 式 的 对 应 的 伴随 方程 是 


—vAw — (Dw)-u+(D*tu)-wt+Vq=ua-u, 在 Q P, 


divw = 0, 在 QD 中 ， 
(19) 
w=0, 在 Tp 上 ， 
o-n+(u-n)w =0, E Pouw E, 
方程 (19) 称 为 方程 (1) 关于 目标 函数 (5) 的 共 斩 方程 。 
同 理 易 得 方程 (1) 关于 目标 泛 函 (6) WHINE 
—vAw —(Dw)-u+(D*tu)-w+Vq=0, #4 Q 中 ， 
divw = 0, 在 Q 中 ， 
(20) 
w= 0, 在 Tp E 
o-n+(u-n)w = ug- u, 在 Tout E. 


4 ”形状 导数 


有 了 共 轿 方程 , 我们 则 可 以 使 用 函数 空间 参数 化 方法 推导 形状 导数 的 表达 式 [591。 设 V(t,z)， 
t > 0, ze 下 "是 区 域 中 的 速度 场 ， 区 域 Q 通 过 速度 场 了 了 得 到 新 的 区 域 0: := TQ) 
& up: = u(Q1), per = pu) 是 方程 (10) 在 变换 后 的 区 域 0: 上 的 解 


EO; Ut, Pt; P, n) 3 0, Vv (9, n) = Vo(e), (21) 
并 且 符 合 目标 函数 使 得 


J(%) = G(R p,n p, Y), (22) 


inf sup 
(gm EVg5 (N) (yy) EVo(Ne) 
其 中 (ut, pe) 满足 


ð 
Z Feit Pei) + 5G (pb N= V (17) € Vol), B. 


G(s; p p, n)õut + GQ wep pnN)6pt =0, Vv(p,n) E Vg(Mz). (24) 
Ut Pt 

FA V(U 是 关于 t+ 的， 下 面 我 们 需要 消去 缚 把 所 有 的 在 区 域 0, 上 和 任意 一 段 边 界 Ft C 
80, 的 积分 转化 到 区 域 0 上 和 边界 PFC ONL. 

由 公式 


sf F(t, x)dax bia = (/ Ei, Z)dz + L F(t, £)V (t); nds) li i (25) 


F) F(t, x)ds les = (人 tt Z)dz + 人 (Fu rz) + HF(t, x) ) V(t) . nds) F , (26) 
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HPF: [0,7] x RN 一 及 是 足够 光滑 的 函数 ,7 是 工 的 外 法 方向 ,， 吾 := divn， 则 当 t = 
ORT, Qo = Q， 对 任意 的 (yp,7) eV5(D)，(W%7Y) € Vo(Q), A 


1 
AG(Q, 9,7, Y, Y) l=0 = [ew — ua)y'ds + 3 人 Ip — ual? Vnds + a(y’, Y) + aly, Y) 
+ f avel) :e(WWVads = b,a) ~ (0,7) 
on 
= f ydiv pVads — B(W", n) — bf, ) 
an 
一 f Ry div WV,ds + clp’; 9, p) + ec(p; p, p) + c(¢; 9, y) 


Doy)p - WVads — (f, 4") — -UVads, 
ae pp WVads — (fy) pre s 
其 中 WV =V(0)-n, Sa! = faoT; "20 = -Va.Y(0)。 更 进一步 的 ， 将 方程 (1) 和 (19) 的 


dt 


解 (u,p; w,q) HA ERE PER (p,n y, y) TU 


AG(Q; u, p; w,q) = A ju — ual?’ Vnds + J 2ve(u) : e(w)Vads 
-f (qdiv u + pdiv wds+ f (Du)u - wVads 
an an 
-f f-wV,ds +a(u, w) — b(u, q’) — b(w, p) + c(u;u, w’) — (f, w") 
an 


+ f (u—ua)u'ds + a(u', w) — b(u’,q) — b(w, p’) + c(u’; u, w) + c(u; u’, w). 
Q 
HR (14) $p =w, n=q, HRI) FY =u, YY =p， 则 上 式 可 以 化 简 为 
dJ(Q;V) = &G(Q; u, p; w, q) = f K (3N; u, p; w, q)Vnds, (27) 
an 
其 中 
K(8Q; u, p; w, q) = Zlu — ual? + 2ve(u) : e(w) ~ qdiv u — pdiv w + (Duju -w — f -w. (28) 


更 进一步 ， 因 只 需要 对 边界 crp 上 进行 优化 ， 则 只 需要 求 得 天 在 直上 的 值 即 可 ， 且 
在 上 divu = 0, divw =0, u=0, W K TRMA 


五 (00i u, p; w, q) |p = zlu — ual? + 2ve(u) : e(w) — f - w, (29) 
易 得 ， 目 标 泛 函 (6) 也 有 类 似 的 结果 
OV) = IGu piwa) = 人 KoniwmwgWas (30) 
an 


其 中 
0 1 
K(0Q; u, p; w, q) |p = By ~ ual? + aHlu — ual? + 2ve(u) : e(w) —f-w. (31) 
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5 ”形状 优化 算法 


现在 我 们 可 以 得 到 形状 优化 的 算法 0%, 蕊 。 首 先 求 出 状态 方程 和 对 偶 方程 的 解 ， 然 后 求 出 目 
标 泛 函 梯度 方向 天 ， 求 解 得 到 下 降 方向 4d 进行 迭代 最 终 得 到 目标 形状 。 因 为 天 仅仅 是 定义 在 边 
界 上 的 ， 所 以 选择 拉 普 拉 斯 算 子 把 天 延 拓 到 区 域 0 上 得 到 玉 ， 而 在 边界 9 LM, KERER 
要 优化 的 边界 下 上 不 为 0。d = (di, d2) € (EL(Q))2， 对 于 任意 的 V c HQ), 4 





va:vwar- (dJ(Q),Vi), += 1,2, (32) 
Q 


且 易 证 4 为 了 的 下 降 方向 。 
则 形状 优化 的 算法 可 以 总 结 为 : 
BRL k= 0， 给 定 初始 的 形状 9 。 
步骤 2 计算 状态 方程 ， 即 计算 其 变 分 形式 (10) 式 ， 得 到 (ut p") Æ |I) 一 J(Q**)| < 
e， 停 止 ， 区 域 f* 为 所 求 区 域 ， 和 否则 转 到 步骤 3。 
PRI WA (ut, p) WARIAT (19) 或 者 (20) 式 ， 得 到 (w*, qg"). 
步骤 4 计算 出 形状 函数 的 梯度 K = K(O0*; ux, pei we, qx)， 用 下 面 方 程 将 其 延 拓 到 区 
RO EER K: 
—-AK*t=0, 72 OF 中 ， 
BPSK. eT! ote (33) 
K* =0, 在 8Q4/I E, 


其 中 工 是 需要 优化 的 边界 。 
步骤 5 求解 下 式 得 到 下 降 方向 下 


-Adk = K*, Æ OF 中 ， 


ode =0, 在 I* E, i= 1,2. (34) 
dt = 0, 在 ONF/TR E, 


步骤 6 给 定 步 长 ak > 0， 形 状 优化 9Q* = OOF! + ardt 得 到 新 的 区 域 Q*, k= 二 上 十 1。 转 
到 步骤 2。 
注 在 这 里 我 们 取 搜 索 步 长 ak 三 a > 0 为 常数 。 


6 数值 算 例 


下 面 我 们 给 出 数值 算 例 验 证 算法 的 正确 性 和 稳定 型 。 取 Reynolds 数 Re = 800, f = 0。 
KN = (0,8) x (0,2)， 图 2 中 实 线 所 示 ， 或 Qa = N/B, BAP WEE (1.5, 0) 半径 为 0.5 的 
上 半圆 ， 图 2 中 虚线 所 示 。 对 于 每 个 算 例 分 别 取 Q A N 作为 目标 形状 来 进行 形状 优化 。 算 
例 1 和 算 例 2 用 来 验证 对 于 两 个 目标 函数 的 算法 的 正确 性 ， 选 择 的 目标 速度 wa 由 在 目标 形状 
下 有 限 元 法 计算 求 得 。 而 算 例 3 要 验证 算法 的 稳定 性 ， 在 目标 速度 wa 上 加 上 了 一 定 的 扰动 。 
如 图 2 虚线 所 示 入 流 边界 Tin 的 边界 条 件 g = (91,92), 9 = y2- y), 92 = 0。 数 值 算 例 使 
用 FreeFem 十 +2.21 软件 进行 计算 图 ， 采 用 P2- _P1 元 。 使 用 计算 机 的 配置 CPU X Pentium 2.4 
GHz， 内 存 为 512MB。 在 实际 工程 计算 中 wu 是 物理 测试 得 到 的 结果 ， 本 文中 取 ua 为 在 精确 区 
域 中 有 限 元 计算 的 结果 。 
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图 2: 灰色 区 域 是 Ql1， 和 白色 的 半圆 型 是 区 域 ， 整 个 算 形 区 域 是 N 


算 例 1 取 目 标 函 数 为 式 (5) 式 。 

图 3 和 图 4 是 以 ?az 为 目标 形状 ， 而 Qi 为 初始 形状 的 优化 结果 。 图 中 显示 的 为 迭代 4 各 步 时 
的 情况 ， 耗 时 843.172 s。 

图 5 和 图 6 以 是 Qi 为 目标 形状 ，Q22 为 初始 形状 的 优化 结果 。 图 中 显示 的 为 迭代 4 各 步 时 的 
情况 ， 耗 时 913.145s。 





一 一 一 ~、 
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图 3: AB: 形状 优化 结果 ; GA: RAR HRESRAPRHAAR 



























































图 4: 水 平方 向 的 速度 wl 图 


算 例 2 取 目 标 函 数 为 (6) xX. 

图 7 和 图 8 是 以 Ql 为 目标 形状 ，Q1 为 初始 形状 的 优化 结果 。 图 中 显示 的 为 迭代 40 步 时 的 
情况 ， 耗 时 867.753 s。 

算 例 3 ”加 扰动 的 算 例 。 

在 实际 计算 中 ， 目 标 速度 函数 wa 常常 由 仪器 测试 得 到 ， 总 不 能 够 十 分 精确 ， 会 有 一 定 的 误 
差 ， 则 要 求 算 法 具有 一 定 的 稳定 性 。 本 文中 为 了 验证 算法 的 稳定 性 ， 讨 论 wa 加 上 3% 的 扰动 之 
后 的 函数 iia 的 形状 优化 情况 ， 则 目标 函数 (5) 变 为 7(Q) = $ fy lu 一 iial?dzx。 
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图 9 和 图 10 是 以 9i 为 精确 区 域 ，92 为 初始 迭代 区 域 迭 代 的 情况 。 图 9 的 左 图 显示 的 为 达 
代 40 步 时 的 情况 ， 耗 时 897.474s。 图 9 的 右 图 是 目标 函数 (9) = 4 folu 一 talde 与 迭代 步 数 
的 关系 图 ( 实 线 和 “十 ”表示 ) ， 还 有 (5) 与 迭代 步 数 的 关系 以 作 比 较 (黑色 点 )。 
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图 5: AB: 形状 优化 结果 ; 右 图 : 和 迭代 目标 函数 值 与 迭代 步 数 的 关系 
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图 7: AB: 形状 优化 结果 ; BA: 和 迭代 目标 函数 值 与 选 代 步 数 的 关系 





图 8: 水 平方 向 的 速度 un 图 
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图 9: AB: 形状 优化 结果 ; GA: 和 迭代 目标 函数 值 与 迭代 步 数 的 关系 














图 10: 水 平方 向 的 速度 wl 图 


7 结论 


算法 求解 混和 边界 条 件 的 不 可 压缩 流 的 形状 优化 问题 是 可 行 的 。 基 本 上 可 以 的 到 期 望 的 形 
状 ， 目标 函数 也 基本 上 收敛 到 0。 

第 一 个 目标 函数 1 即 (5) 比 第 二 个 目标 函数 即 (6) 收敛 速度 快 而 且 收敛 效果 明显 ， 得 到 的 区 
RO 也 更 加 接近 精确 解 。 因 为 (5) 中 求解 的 速度 与 目标 速度 va 在 整个 区 域 上 进行 比较 ， 能 得 到 
更 多 关于 目标 速度 wa 的 信息 ， 而 (6) 仅仅 是 在 边界 上 进行 比较 。 

该 算法 具有 一 定 的 稳定 性 ， 当 目标 速度 函数 wa 加 上 一 定 的 扰动 后 仍然 可 以 得 到 比较 有 效 的 
结果 。 

更 进一步 这 些 算法 可 以 应 用 到 非 定常 的 Navier-Stokes 方程 ， 而 且 也 可 以 应 用 到 工业 中 去 。 
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Function Space Parameterization Technique for Shape Optimization 
Problem of Incompressible Fluid 
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Abstract: The shape optimal design for fluid is very important in applications. For a fluid in a region 
described by a Navier-Stokes equation, our cost fuction is a functional of states of the fluid. Based on 
the adjoint method and a parametrization method, a formulae of the shape derivatives is established 
for the problem. Then a conjugate-graduate algorithm is proposed. The numerical examples show the 
effectiveness and stability of the algorithm. 
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